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Abstract 

How can repulsive and attractive forces, acting on a conservative system, create stable traveling 
patterns or branching instabilities? We have proposed to study this question in the framework of 
the hyperbolic Keller-Segel system with logistic sensitivity. This is a model system motivated by 
experiments on cell communities auto-organization, a field which is also called socio-biology. We 
continue earlier modeling work, where we have shown numerically that branching patterns arise for 
this system and we have analyzed this instability by formal asymptotics for small diffusivity of the 
chemo-repellent. 

Here we are interested in the more general situation, where the diffusivities of both the chemo- 
attractant and the chemo-repellent are positive. To do so, we develop an appropriate functional 
analysis framework. We apply our method to two cases. Firstly we analyze steady states. Secondly 
we analyze traveling waves when neglecting the degradation coefficient of the chemo-repellent; the 
unique wave speed appears through a singularity cancelation which is the main theoretical difficulty 
This shows that in different situations the cell density takes the shape of a plateau. 

The existence of steady states and traveling plateaus are a symptom of how rich the system is 
and why branching instabilities can occur. Numerical tests show that large plateaus may split into 
smaller ones, which remain stable. 

Key-words. Keller-Segel system, hyperbolic system, traveling waves, branching instability, cell com- 
munities. 
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1 Introduction 



We study a form of the hyperbolic Keller-Segel system with logistic sensitivity for a chemoattractant 
and constant sensitivity for a chemo-repellent, given by the set of equations 



dtn + div 
-D c Ac + 



fjL c n 1 



Vc — fis n ^S 



0, 



a c n . 



(1) 



[ d t S-D s AS + ^ = a s n. 
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This system is reminiscent of a large class of models used with biomedical motivations to represent 
the auto-organization of cells that are able to produce attractive and repulsive chemicals [H [18] , but 
also in other areas of population dynamics [H [TS]- The logistic sensitivity /x c (l — n/n max ) takes 
into account a volume filling (or quorum sensing) effect, i.e., a reduction of the cell response to the 
chemo-attractant (whose concentration is denoted by c), which prevents overcrowding [1(J|. \T§\ I2(J| . 
The special form at hand has been proposed in [5j as a reduced model for a more detailed system to 
study complex patterns as the dendritic ramification of Bacillus subtilis, recently obtained with high 
nutrient experiments in [111 [T2"| [Tl] , whereas pattern formation based on local nutrient depletion is 
also possible [TU El C£7]. It includes a chemo-repellent of concentration S (that can be interpreted 
as the effect of surfactin) with a constant sensitivity us- The model neglects diffusion (i.e. random 
motion) of the cells. This has been proved to be mathematically correct, due to the logistic sensitivity, 
in [6] (see also [2] for an earlier work) . Reaction-diffusion models are used for the chemicals with 
a quasi-stationarity assumption for the chemo-attractant. The chemicals diffuse with diffusivities D c 
and Ds, they are degraded with relaxation times r c and t$, and they are produced by the cells with 
rates a c and as- 

From [U [2j [7] we know that, when 5 = 0, the nonlinear term n(l — n/n max ) causes sharp fronts 
that connect alternatively the states n = and n = n max (see [9] for this terminology). The repellent 
force V5 can generate surprising dynamics of the plateaus and branching instabilities may occur. 
This was shown in [SJ, and the instability could be analyzed for Ds small, because the limiting system 
with Ds = can be recast as an hyperbolic system according to a method introduced in [13] . Then, 
stability/instability of discontinuities can be seen as a transition from shock to rarefaction waves. 

In the present paper, we are interested in the existence and branching instabilities of traveling plateau 
solutions of ([I]) for general diffusion coefficients. These are naturally obtained in one dimension and 
therefore we focus on the nondimensionalized system 



For the nondimensionalization, ts has been chosen as the reference time and n max as the reference cell 
density. The reference values for the chemical densities have been chosen such that the scaled versions 
of the sensitivities fi c and fis are equal to 1. All quantities in are dimensionless. This is also true 
for the parameters D c , Ds, a c , and as, where we point out that the scaled diffusivities D c and Ds 
actually represent (in unsealed notation) D c t c /1 2 and Dsts/1 2 , respectively, where I is the reference 
length. 

We develop a functional analytic framework for the study of traveling waves. Since it is simpler to 
introduce it for steady states of ([2]), we address this issue in Section [2l We prove the existence of a 
family of steady states characterized by the size of the plateau, when it is small enough, or when Ds 
is close enough to D c (Section [3]). Numerical tests illustrate that the smallness condition is necessary 
for stability. In Section SJ we show that steady states are replaced by traveling plateaus, when 
degradation of the chemo-repellent is neglected (i.e. ts = oo in ([I])). The method is an extension 
of the functional analytic framework for steady states, where the propagation speed is determined 
naturally by a singularity analysis. The total number of cells (or the size of the plateau) defines a 
family of traveling waves with different speeds. 




< 



—D c d 2 c + c = a c n , 
k d t S-D s d 2 x S + S = a s n. 



(2) 
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2 Steady states with small total mass 



The existence of stationary states for system ([2]) is both the most natural question and the easiest to 
illustrate the method we use throughout the paper. Therefore we begin with this issue. We first state 
a theoretical result that involves a smallness condition. Then, we confirm with numerical results that 
this smallness condition is necessary. 
We shall look for solutions of the steady state system 

d x (n(l — n)d x c — nd x S~) = , 

—D c d x c + c = a c n , (3) 
-D s dlS + S = a s n, 

complemented with the boundary conditions 

ra(±oo) = c(±oo) = s(±oo) = . (4) 

Integration of the first equation then gives vanishing flux: 

n[(l-n)d x c-d x S] =0. 

We shall be interested in plateaus of the cell density with sharp boundaries, such that n jumps between 
n = and 

- = l-g>0 (5, 

(satisfying the Rankine-Hugoniot jump conditions). The convexity of the flux function is determined 
by the sign of d x c. The above jumps satisfy the entropy condition if either n = on the left, n > 
on the right of the jump, and d x c > at the jump point; or n > on the left, n = on the right, and 
d x c < 0. 

Theorem 2.1 Assume that 7 := < 1. Then, for L small enough, there exists a unique entropy 

solution of d3J) of the form 



n(x) 



■ for x^ [0,L], 

, l- 7 + 0(L) for xe(0,L), 

with n £ C(0,L) and symmetric around x = L/2. 



For the existence of steady state plateau solutions, we need the effect of the chemo-attractant to 
dominate the effect of the chemo-repellent. In terms of the original unsealed parameters, 7 = a ^ s 
holds. Thus, the strength of the attractive (respectively repulsive) effect is measured by the product 
of the production rate and the sensitivity divided by the diffusivity of the chemical. 

One can understand the occurrence of the free parameter L as the result of mass conservation in 
the dynamics ([2]). Supposedly there is a one-to-one relation between L and the total number of cells. 
This is also the way plateaus of different size are produced in the numerical examples below. 
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Proof of Theorem 12.11 

The difficulty in using (0) for the computation of the cell density is to control the points where d x c 
vanishes. By the symmetry assumption, one such point is at x = L/2, and we shall prove that for n 
of the form given in the theorem, it is the only one. 

It will be convenient to rescale the problem by x — > Lx. Then, the support of n is given by [0, 1], 
where ([5]) still holds, and 

( -^d 2 c + c = a c n, 
\ -^rd 2 x S+S = a s n. 
With the boundary conditions explicit representations of d x c and d x S can be computed for 

d x c = -^(J\ L ^n{y)dy- j\ L ^n(y)dy^, (7) 

d x S = -^^^j\ L s(y--) n ( y )dy- j\ L s(--y)n(y)dy^ , (8) 
with L$ = L/^Ds, L c = L/^/D c . Thus the formula ([5]) can be seen as the fixed point equation 

with 

F[Li, n](x) = - [ e Ldy - x) n{y)dy + [ e L ' (x - y) n(y)dy , < x < 1 , i = S, c. (10) 

Jo Jx 

Using the symmetry at x = 1/2, i.e., n(x) = n(l — x) for x S (0, 1), we restrict to x G (0, 1/2) and 
write 

fi/2 ex 
F[L h n]{x) = e L * x J (e~ LiV + eT^M) n(y)dy + (e^ 1 ^ - e~ Ux } J e Uy n(y)dy. (11) 

As the next step, the limit as Lj — > is taken out, and the zero at x = 1/2 is eliminated by 

F[Li,n](x) 



1/2 - x 

with 



2n(x) + Lif[Li,n](x) , 



1/2 

I /-1/2 e Li(x-y) _|_ e Li(x+y-l) _ 2 g — x) _ e ~LiX rx 



n(x) = - J niuultj. 



f[L ^ n]{x) = Tj2—~xj x ~ ^ Mv)dv + - Li(1/2 _; } ^ e^n(y)dy. 

Since L = implies L, L = 0, the unique solution of ([9]) in this case is given by 

n(x) = no := 1 — 7 . 
We shall apply the Banach fixed point theorem in the ball 

B:=\n€ C([0, 1/2]) : ||n - n ||oo < , (12) 



of the Banach space C([0, 1/2]) (equipped with the supremum norm || • ||oo)* The essential observation 
is collected in the following result. 
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Lemma 2.2 For < Lj < L, the linear mapping n i— >■ /[Lj,n] on C([0, 1/2]) is bounded with a bound 
only depending on L. 

Proof. It is easily seen that both 

e Li(x-y) _|_ e Li(x+y-l) _ 2 g — a;) _ g— LjZ 

and — — ; r 

U Li{\/2-x) 

are uniformly bounded for < x, y < 1/2, < Li < L with a bound C(L), implying 
□ 

Corollary 2.3 For n G B and for L small enough (recalling L c = L/\JTT C , L$ = Lj VAs ), 

2n + LJ[L c ,n] > — ^ 

and 

irM , . L c f[L c ,n] - L s f[L s ,n] 

n H- T[n\ = 1 - 7 + 7 „_ - T , rr ; 

2n + L c /[L c ,nJ 

maps -B into itself and is a contraction with respect to || • ||oo with a Lipschitz constant proportional to 
L. 

This concludes the proof of Theorem 12.11 

Formal asymptotic expansion — shape of the plateau 

In this section the first few terms in an asymptotic expansion of the solution of ([9]) will be computed. 
This will shed light on the shape of the non-constant correction of the cell density plateau. Some of 
the necessary computations are rather lengthy and will only be outlined. We start with the Taylor 
expansion F[Li,n] = F [n] + LiF^n] + L 2 F 2 [n] + O(Lf) with 

F 3 [n](x) = - f X ^LZpL n{y)dy+ f ^yL n{y)d y . 

Jo r- Jx r- 



This in turn leads to 



F[n] = T [n] + LT x [n] + L 2 F 2 [n] + 0(L 3 ) 



with 



T [n] = 1-7, Ji[n] =7 



fD e VDsJ F [n] 
1 1\ F 2 [n] 7(1 I \ Fi[n\ 2 



D c D S J F [n] ^D~ c \^ yfJTJ F [n} 2 ' 
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Substitution of the ansatz n = no+ Ln\+ L 2 n 2 +0(L 3 ) into ([9]), re-expansion, and equating coefficients 
of powers of L then leads to 




n = 1 - 7 , 

7 

m = 7"i [no] = - 



n 2 = DJi[no]ni + 7 2 [n ] = n 2 + | f ^- - -^J (s - 1/2) 2 , 

where DT\ denotes the Frechet derivative of T\ and n 2 is a (explicitly computable) constant. The 
0(L)-correction term n\ is constant. It is negative for D c < D$ and positive for D c > D$- The first 
non-constant correction occurs at 0(L 2 ). It is convex for D c < D$ and concave for D c > D$- This 
agrees qualitatively with the numerical results of Section [3l 

Finally, we mention that it is a standard procedure to extend our rigorous results in order to justify 
the asymptotic expansion for n in the sense that the error 0(L 3 ) can be estimated in C([0, 1/2]) by 
CL 3 . 



Numerical experiments 

We carried out numerical tests that illustrate the analytical results and indicate that large plateaus 
may be unstable, depending on the relation between D c and D$- 

We obtained numerical steady states as the limit for large times of a modified dynamics (where 
also the chemo-repellent is determined from a quasistationary problem), which we have chosen for its 
simplicity: 

dtn + d x (n(l — n)d x c — nd x S) = , 
< -D c d^c + c = o c n , (13) 
-D s dlS + S = a s n. 

We discretize the hyperbolic equation for n by the Enquist-Osher finite volume method, which is 
conservative and can capture the shocks on both sides of the plateaus (see [1] for a recent introduction 
to the subject). The elliptic equations for c and S are solved by a finite difference method. We indeed 
obtained that after some transient the numerical solutions converge to a steady state. 

We present three families of results in Figures [H [2] and [3j In these pictures, the solid and dashed 
lines on the sub-figures at the top represent n and c respectively, while the bottom sub-figures depict 
S. 

First we illustrate Theorem 12.11 (that is L small) in Figure HJ Here the computational domain is 
[0, 6] and the initial density is an almost centered indicator function: 

n° = 1[ 3 ,3.5] • 

We observe the different shapes of the steady states supporting the formal asymptotics above. When 
D c > Ds, the positive part of n is concave, when D c = D$ it is flat, and when D c < D$ it is convex. 
The leading order approximation for small L suggests the relation M ~ L(l — 7) between the total 
mass M and the width of the plateau. Since M and L are not very small in these computations, this 
approximative relation should, however, be corrected by higher order terms. In particular, note that 
in the simulations leading to the right picture, 7=1 holds, such that the leading order term does not 
provide any contribution to the total mass. 
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Figure 1: The steady state solution of ([3]) with total mass M = 0.5 and different chemical production 
rates and diffusivities. Left: a c = I, as = 1,D C = 0.1, D$ = 1; middle: a c = l,as = 0.1, D c = 
0.1, D s = 0.1; right: a c = l,a s = 0.1, £> c = 0.1, L> 5 = 0.01. 



Secondly, we test for L big. The numerical results for a c = 1, as = 1, D c = 0.1, Ds = 1 are presented 
in Figure [21 with initial data corresponding to M = 1, namely 

n°(x) = I[2.5,3.5] • 

It seems that the total mass is too large for a one-plateau steady state to exist. The initial plateau 
splits into two smaller ones which appear to be stable. 

It appears numerically that, when 7 < 1, D c > Ds, we always reach a steady state solution with 
concave cell density in the plateau, no matter how large the total mass is. As an illustration we show 
numerical results with M = 2 and the initial density 

n°(x) = l [2i 4] • 

The results are shown in Figure El 

3 Steady states for almost equal diffusion lengths 

Obviously, for D c = Ds, i.e. L c = Lg, and for 7 = as/a c < 1 the steady state problem Q has the 
constant solution n = uq = 1 — 7 (for arbitrary L > 0). Note that, in terms of the original unsealed 
parameters (as occurring in ([T])) the equality of the scaled diffusivities means equality of the quantities 
\JD c t c and \fD~s~rs, which can be interpreted as diffusion lengths, i.e. the average distance a molecule 
diffuses, before it gets degraded. 

When considering the dynamics (fT3|) with D c = Ds, it is obvious that cas = Sa c holds. Therefore 
(fT3l) is equivalent to 

d t n + d x (n(l - 7 - n)d x c) = , 
—D c d x c + c = a c n . 

For 7 < 1 this is exactly the problem analyzed in [7\, where strong arguments for the stability of 
one-plateau solutions are given. 

In the remainder of this section we carry out a perturbative analysis to prove existence of a plateau 
steady state for small values of Lg — L c , and we illustrate the qualitative behaviour of the solution by 
a formal asymptotic expansion in terms of the perturbation parameter. 
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Figure 2: The instability for D c < D$ when M = 1. We observe that the plateau splits in two pieces, 
each of which stabilizes after they are well separated. In these four pairs of figures, the solid and 
dashed lines on the top sub-figure represent n and c, respectively, and the bottom sub-figure shows S. 
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Figure 3: The stability of the plateau when D c > D$, even with M = 2. The parameters are 
a c = 1, as = 0.1, D c = 0.1, Dg = 0.01. The solid and dashed lines on the top sub-figure represent n 
and c respectively and the bottom sub-figure shows S. 
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Existence of a steady state 

Our approach is completely analogous to the preceding section. We fix L c > 0, introduce the pertur- 
bation parameter 5 := L$ — L c , and use the decomposition (fTTjk 

F[L c + 5,n] 



1/2 -x 

with G[n] = F[L c ,n]/(l/2 - x) and 



G[n) + 5g[S,n], 



1 f 1 / 2 

9[^n\ = -— / ^(8,x,y)n(y)dy + I tp(8,x,y)n(y)dy , 



1/2 -x 
where 

e (L c +S)(x-y) _|_ e (L c +6)(x+y-l) _ e L c (x-y) _ e L c (x+y-l) 

<p(8,v,y) = 7 > 

e (L c +8)(x+y-l) _ e {L c +S)(y-x) _ e L c (x+y-l) + e L c {y-x) 

mX > y) = Wj2^x) • 

With this notation, the fixed point problem ([9]) reads 

n = !-,-*ff (15) 

So we need uniform (in 5) boundedness of the linear maps g and G, as well as boundedness away from 
zero of G[n]. 

Lemma 3.1 Let 7 = a s /a c < 1 and let 

B:={n£ C([0, 1/2]) : ||n - 1 + 7IU < C B ) , with C B = ^ ~ l)e ^ . 

2(2 — e ~ ljc / z ) 

Then 



G[n] > 



(l- 7 )e- Lc / 2 (l 



e 



2L C 

for n G -B cmc? i/ie map G is bounded with respect to || • ||oo* 

Proof. The boundedness follows immediately from the representation 

p L c x pl/2 . . p -L c x _ p -L c (l-x) px 

G[n]{x) = 1/2^ J V + e ~ Lcil ~ y) ) n (y) d V yj—^ J e L ^n(y)dy . 

Since Cb < (1 — 7)/2, n > for n G B. Therefore, for x G (0, 1/2), G[n] is the difference of two 
positive terms, and we can estimate 

p L c x(-\ _ ~, _ n n \ /"1/2 , s -L c x _ —L c (l-x) rx 

e (1 7 Ob) I f„-L c y, p -L c (i-y)\i_e e / L 



G[n](x)> J (e-^ + e-^ L -y))dy — (l- 7 + C B ) J e^dy. 

Evaluation of the integrals leads to 

0| " ] W £ ft^d-l)) [1 - 7 - Ob - (1 - 7 + ft) (1 - a—)] , 

where the Bernoulli function Be(z) = z/(e z — 1) is positive and strictly decreasing. By setting x = 
in its argument and x = 1/2 in the bracket, the lower bound from the statement of the lemma is 
achieved. □ 
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Lemma 3.2 For < 5 < 5, the linear mapping n \— > g[S, n] on C([0, 1/2]) is bounded with a bound 
only depending on 5 and L c . 

Proof. The result is a straightforward consequence of the observation that the functions tp and ip 
are bounded for 5 G (0,5], x, y G [0, 1/2). This again follows from the facts that their denumerators 
<p(S, x, y) and, respectively, ip(5, x, y) are smooth functions of their arguments satisfying <p(0, x, y) = 

$(0,x,y)=$(8,l/2,v)=0- □ 

The last two results immediately imply the contraction property of the right hand side of (|15p acting 
on B for small enough 5, which proves the following existence result. 

Theorem 3.3 Assume that j := < 1 and define S = L{D S — D c ). For \5\ small enough 
(with L and D c fixed), there exists a unique entropy solution of ([3]) of the form 



n(x) 



■ for xi [Q,L], 

, 1-7 + 0(5) for x £ (0,L), 
with n G C(0, L) anc? symmetric around x = L/2. 



Formal asymptotic expansion - shape of the plateau 

Similarly to the preceding section we start with the expansion 

F[L C + 5,n]= F[L C , n] + 5F l [L c , n] + 0(5 2 ) , 

with ^ 

A [L c , n] = f e L "^ (x - y)n(y)dy + [ e L ^ (x - y)n(y)dy . 

JO Jx 

This leads to the asymptotic expansion for the cell density: 

„ = i- 7 - ^' 1 --;W ). 

F[L C , 1 - 7] 

A straightforward computation gives 

Differentiation and the inequality e n — e~ v — rj > for rj > imply that the function in the parentheses 
is strictly increasing for x G (0, 1/2) and strictly decreasing for x G (1/2, 1). Therefore n has the same 
property for 5 > 0, i.e. D$ < D c , and the opposite for D$ > D c . These are the same qualitative 
results as in the preceding section. 



4 Existence of short traveling plateaus 

So far, we have proved the existence of stationary solutions, which can be viewed as traveling waves 
with zero velocity. In this section, we use a similar route to establish the existence of some non-zero 
velocity traveling plateaus. These exist when degradation of the chemo-repellent is neglected, i.e. 
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Ts = oo in dU). With an appropriate nondimensionalization and a reduction to one dimension, the 
system becomes 

d t n + d x [(l - n)nd x c - nd x S] = , 

—D c d x c + c = a c n , (16) 
^ d t S - D s d 2 x S = a s n. 

We consider plateaus of length L and with speed a. They are defined as functions of the traveling 
wave variable x — at, which for simplicity is again denoted by x: 

—ad x n + d x [{\ — n)nd x c — nd x S] = , 
-D c d x c + c = a c n , 
-ad x S — Dsd x S = asn . 
Again, we restrict our attention to solutions satisfying 

n{x) > , x G (0, L) , n = , else. 
We rescale space as x — > Lx and the wave speed as a — > La and obtain for x 6 M 

—al?d x n + d x [(l — n)nd x c — nd x s] = , 



(17) 



-jfid x c + c = a c n , 
-ad x S - ^d 2 x S = a s n , 



with 



(18) 



(19) 



n > in (0, 1) , n = else. 

System (fT8|) - (fT9]) is defined in the whole space x E R and we complete it with boundary conditions 
for c and S: 

c(±oo) = 0, S(-oo) = , S"(+oo) = 0. (20) 

The wave speed and the far field value Sqq of the chemo-repellent are considered as unknown and part 
of the solution. 

Theorem 4.1 For 7 = ° S p° < 1 and L small enough, there is a unique solution of (fT8|) - ([20]) . such 
that n 6 W l '°°(Q, 1) and 

n = l- 7 + 0(L) in (0,1), <r = ~ 7) +Q(£), = 2Z> S + O(L) . 

Furthermore, c is concave in (0, 1) and S is non-increasing on R. 

The method of proof extends, with additional technicalities, that of Section [2j After integrating the 
equation for n, we find the following formulas for the solutions of (|18|) — fj2Q p : 



n 



1 



ah? + d x S 



Or. 



x G [0,1], 



d r c = 



oLc(y~ x ^n(y)dy+ / e L ^ x -^n{y)dy 



x e R, 



d x S 



-a s L 2 s T e° L2 s^n{y)dy, xGR, 
Jo 



(21) 
(22) 
(23) 
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where, again, the notation L c = L/^fD~ c , L$ = L/^/Dg has been used. If the wave speed a was 
known, the right hand side of (l2Tj) could be (after substitution of ([22]) . (j23j) ) considered as a fixed 
point operator for the computation of n. The difficulty is that a is not known a priori. The new 
ingredient compared to the previous sections is that a is used to create a zero of the denumerator at 
the same position where it occurs in the denominator. This principle to find a traveling wave speed 
seems to be new. 

Formally, the procedure (whose feasibility will have to be proven) is as follows: Given n, find 
= such that F[L c ,n](xo) = 0, with F[L c ,n] as defined in (fTU]) . Then determine a = a[n] 

such that 

a = ^f f ° e aL2 s iy - xo) n(y)dy . (24) 
D s Jo 

Now the cell density in the plateau can be determined as a fixed point of 

2Dc a[n]-?jLfi<rm<v-*)n(y)dy 

™ x) = 1 - — ■ (25) 

The next step is to show that, for L small enough, J-[n] is well defined on 

1-7' 



B = in G C([0,1]) : ||n - l + 7 ||oo < 
i.e. for (1 - j)/2 =: n < n < n := 3(1 - j)/2. 

Lemma 4.2 For L c small enough and n G B there exists a unique xq = xo[n] G (0,1) satisfying 
F[L c ,n](xo) = 0. Its dependence on n is Lipschitz: 

\x [ni] - x [n 2 ]\ < C\\ni - n 2 \\oo , 

for all ni,n 2 G B, with C independent from L c —> 0. 

Proof. Existence follows from continuity of F[L c ,n] and from F[L c ,n](0) > and F[L c ,n](l) < 0. 
Introducing Xj = xo[nj], j = 1,2, the equation F[L c ,rij](xj) = can be written as 

[ Xj e L ^ nj (y)dy = e 2L ^ f e^n^dy , j = 1,2. 

JO Jxj 

The difference of these two equations is written in the form 

[e L -y + e Lc{2x 2 -y)\ ni (y) dy + ( e 2L c x 2 _ e 2L c x^ j ^-Uy^ ( y ) dy 

J X2 J X\ 

= r e L ^(n 2 (y)-n 1 (y))dy + e 2L ^ f e~ L ^(n^y) - n 2 (y))dy . 

JO Jx 2 

Now the modulus of the left hand side is estimated from below in terms of \x\ — x 2 \ and the right 
hand side from above: 

(n — 2L c e Lc n) \x\ — x 2 \ < 2e Lc \\m — n2||oo • 

Since the coefficient on the left hand side can be made positive by choosing L c small enough, this 
concludes the proof. □ 
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Lemma 4.3 Let the assumptions of Lemma be satisfied. Then there exists a unique a = o~[n] 
satisfying (|24p with xq = xo[n}. Its dependence on n is Lipschitz: 

\o-[m] - a[n 2 ]\ < C\\ni - n 2 \\oo , 

for all n\,ri2 £ B, with C independent from L$, L c — > 0. Furthermore, 

"SO- ~ 7) e - Ls * =; a < a < - ;= 3a5(l ~ 7) 
2D S -~ ~ 2D S 

Proof. Existence and uniqueness of a positive solution follow immediately from the fact that the right 
hand side of (|24p is positive, bounded, and nonincreasing as a function of a. The upper bound is a 
consequence of n £ B, and so is the lower bound (in the derivation of which the upper bound is also 
used). Similarly to the proof of Lemma 14.21 for ni,n 2 G B the difference between the corresponding 
cr-equations can be written as 

pX\ f*X\ 

o-i - 0-2 - ^- I ( e ^ L s(y-^i) - e ^L s {y-x^ ni ( y ) dy = I e^ Ls ^- x ^m{y)dy 

X2 



DsJo V J ^' D s 



+TT / {e a2Ls{v - Xl) - e^ Ls{y - x A ni {y)dy + / e^ Lsiy ~ X2) ( ni (y) - n 2 (y))dy . 
L>s Jo v ' B> s Jo 

It is now straightforward to estimate 

\a x -a 2 \<C {\x x - x 2 \ + \\rii - n^^) , 

and to complete the proof by using Lemma 14.21 □ 
With the definitions of xo[n] and cr[n], the fixed point operator can be written as 

F[n] = g[x [n],a[n],n] 

with notations where the singularity cancelation at xq appears more clearly 

r\ l 1 o I \ x 0,n]+ g[xQ,<r,n) f x ° 

g[x ,a,n] = l-2 7 — — — -, I[x ,n](x)= n(y)dy , 

2L[x ,n\ + f[x ,n\ J x 

xo / \ r x 



g[x ,a,n](x) = J {e vL ^~^ - 1^ n{y)dy - J - l) n(y)dy 

f'XO . . f'X . . 

= J (e^y-**) -Vjn{y)dy + J ^e aL s^- x ^ - e^y-^j n(y)dy , 

/ x ° / \ 
( e L c (y-x ) + e L c (x -y) _ 2 j n ^ dy 

^L c (y-x ) _ e L c (y-x)^ n ^ dy + / ^L c (x~y) _ e L c (x -y)^ 



Since / and g vanish for L = (=>• L$ = L c = 0), the constant no = 1 — 7 is the only fixed point in 
this case. The following properties of /, / and g are obtained by straightforward computations. 

Lemma 4.4 For n £ B there exists a constant C independent from L — >• 0, such that 

1 L[x ,n](x) 3 

2 (1 " 7) ^ x -x ^2 (1 ~ 7) ' (26) 
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\g[x ,a,n}(x) - (x - x )g[x ,a,n}'(x)\ < CL 2 (x - xq) 2 , (27) 
\f[x ,n](x) - (x - x )f[x ,n]'(x)\ <CL(x-x ) 2 , (28) 
\f[x ,n]'(x)\<CL, \g[x ,a,n}'(x)\<CL 2 . (29) 
Moreover C only depends on Halloo — and not on n f . 

Since we already have the Lemmas 14.21 and 14.31 we only need to examine the dependence of Q[xq, a, n] 
on its arguments for proving the contraction property of J-. Unfortunately, it will turn out that the 
Lipschitz constant of Q as a function of xq involves the derivative of n with respect to x. Therefore, 
we shall need a stricter definition of the set, where the fixed point iteration is carried out. 

Lemma 4.5 There exists a positive constant k such that, for L small enough, the fixed point operator 
T maps both B into itself and the set 

B = {n€B: n € W 1|DO ([0, 1]), < kL] 

into itself. 

Proof. We use the alternative representation 

x-r i , 2g[x Q [n],(r[n},n} - f[x [n],n] , . 

^ N = 1 - 7 - 7 2/[*oN,n] + / M n],n] ^ 

and the controls 

I[x ,n](x) > lzl 2J(g) + f(x) > 1 - 7 CL (31) 

Xq — X ~ 2 ' Xq — X ~ 2 

As a first step, our previous result, together with the estimate for a[n] in Lemma 14.31 after cancel- 
lation of x — xq in the denominator and the denumerator, gives 

CL 



\F[n] - 1 + 7 | < 



1-j-CL' 



This implies that, for small enough L, T maps B into itself. The second step is to compute the 
x-derivative of the fixed point operator: 

r[ y 9 I'{2g-f)-I{2g'-f')+f'g-g'f 

g[xQ,a,n\ = 2 7 (21 + f) 2 ^ ' 



For estimating this term, we need (|27p . (|28[) and the corresponding property of /. This is a first time 
when the derivative of n enters: 



\l[xo,n](x) - (x - x )I[xo,n]'(x)\ < 



n(y)dy — (xq — x)n(x) 



< ~(x - x ) 2 \\n'\ 



With these properties, for n £ B, the modulus of the denumerator in (|32p can be estimated from 
above by CL(1 + kL)(x — xq) 2 . On the other hand, the denominator can be estimated from below by 
(1 — 7 — CL) 2 (x — xq) 2 . Thus, the bound on the derivative is preserved by J 7 , if 

C(1 + kL) 

< K . 



(l- 7 -CX)2 

This holds for any k > C(l — "f)~ 2 and small enough L. □ 
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Lemma 4.6 For L small enough, < xoi, X02 < 1> o"i? o"2 satisfying the bounds in Lemma \4-3[ and 
rix, n 2 G -B, 

||C7[a;oi,<7i,ni] ~ £[^02, o" 2 , n 2 ] < CL (|x i - X02I + l a i ~ °2| + ||ni - n 2 \\ 00 ) , 
wrai/i C independent from L. Moreover C only depends 

Proof. For analyzing the dependence of Q[xQ,a,n] on xq, it is convenient to observe the iden- 
tity G[xo, a, n](x) = G[x,a,n](xo), following from the skew symmetry of I[xo,n](x), f[xo,n](x), and 
g[xo,a,n](x) with respect to x and xq. It implies 

- — Q[x ,o-,n](x) = Q[x,a,n]'(xo) . 
dxQ 

A bound of the form CL of this quantity has been shown in the proof of the previous lemma. 
From the definition of g[xo, o~, n] it is obvious that the derivative with respect to a is 0(L 2 ). 
Finally, I[xo, f[xo, -](x) and g[xo,a, -](x) are linear functionals with the obvious bounds 

00 

Considering again the representation (|30jl and the lower bound (I31|) . the proof is completed. □ 

Combining the Lemmas l4.2l 14.31 14.51 !4.6l shows that, for L small enough, T : B — > B is a contraction, 
completing the proof of the existence and uniqueness statement of Theorem 14.11 The limit as L — > 
of n follows from the form (|30p of the fixed point operator. Using this limit in the limit of the equation 
F[L c ,n](xo) = shows that the limit of xq is 1/2. Then the limit of a is obtained from L$ — > in 
AMD ; and the limit of the far- field value Soo is derived by using ([23]) : 

POO 

Soo = S(-oo) = 5(0) = - / d x S{x)dx -> 2D S . 

•/ 

This completes the proof of Theorem 14.11 



I[x ,n] 



X — Xq 



< I nil 



/[x ,n] 



x — x 



< CLHnll 



# |a;o,cr, n] 



X — Xq 



5 Numerical simulation of traveling plateaus 
Algorithm 

In the previous section, the existence of traveling wave solutions of system (|16p in the form of short 
enough cell density plateaus has been proven. In the following, numerical simulations of system (jl6[) 
will be presented, indicating the necessity of the shortness assumption for the dynamic stability of 
traveling plateaus. 

In order to obtain fast enough convergence to a traveling wave, while it moves through the finite 
computational domain [0, A], an approximation of the wave, corresponding to the limit L — > 0, is used 
as initial condition. 

More precisely, the initial cell density n(t = 0) is chosen as 

n = (1 -7)l(i,i+L) , (33) 
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where we recall that 7 = ■ With this cell density, the maximum of the chemo-attractant concen- 
tration occurs at xq = 1 + L/2, and the velocity a is obtained by solving the nonlinear equation (|24p . 
Next, an initial datum Sq for the chemo-repellent concentration is computed by solving the equation 
for 5 in ([HI) (see also 

Since the numerical scheme is restricted to a finite computational domain [0,A], we need to define 
appropriate boundary conditions. Simulation times T > are chosen such that the set {x, n(t, x) > 0} 
stays away from the boundary points x = 0, A for < t < T and, since the equation for n is hyperbolic, 
it is enough to use zero entering flux boundary conditions. 

For the chemical concentrations, Robin boundary conditions are used, which are satisfied exactly 
by traveling plateau solutions: 

sfW c d x c(t, 0) = c{t, 0) , tJWc 8 x c{t, A) = -c(t,A), (34) 

d x S(t,0) = 0, D s d x S{t,A) + aS{t,A) = 0. (35) 

The value of a needs to be updated for each time step. This requires the computation of a[n] as 
described in the previous section, where the cell density n from the previous time step is used. The 
nonlinear equations for xo[n] and a[n] are solved by the Newton method. 

As in Section [2j the equation for n is discretized by the Enquist-Osher scheme and the quasi- 
stationary equation for c by finite differences. For the parabolic equation for S a finite difference 
space discretization is used with time implicit treatment of the diffusion term. 

Numerical results 

Effect of L: In all our simulations, the computational domain is [0,6], i.e. A = 6. For the fixed 
choice of parameters 

a c = 1, a s = 1, D c = 0.1, D s = 1, 

we take different values for the length of the initial plateau. For L = 0.1, 0.4, and 0.7, the evolution of 
re, c, and S is plotted in Figures [U and[H respectively. We observe that, when L = 0.1 or 0.4, the 
solution converges to a plateau traveling wave. Note the differences in the plotted times between the 
two figures, indicating the different wave speeds. For the largest value of L, the initial plateau splits 
into two pieces that travel independently (with speeds dictated by their post-splitting lengths). An 
interesting question (we do not have an answer to) is, whether a dynamically unstable traveling wave 
still exists. 

Plateau shapes for different physical constants: Different parameters can give different shapes 
for the plateau. Figure [7] depicts the detailed shape of the plateau depending on the diffusion coeffi- 
cient of the chemical. It seems that when D c (and therefore 7) increases, the plateau becomes larger 
and layers appear on the edges of the plateau. In this case, attraction forces decrease and therefore 
cells diffuse more in the middle. 

Convergence of the scheme: We recall that the wave velocity a[n] is computed in each time step. 
The numerical values are displayed in Figure [8] for two different values of L. It seems that the value 
of a computed from the initial data is already very close to its steady state value. Though there are 
small oscillations, when the mesh is refined, the amplitude of the oscillations is reduced, indicating 
convergence of the numerical method. It is easy to see that, regardless of the oscillations, a is almost 
constant. Specifically, its value is close to ^-L(l — 7) in accordance with Theorem 14. 1[ 
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Figure 4: The solution of ()16|) with L = 0.1 for three different times. Top: The solid and dashed lines 
are n and c, respectively. Bottom: the evolution of S. Left: t = 0, middle: t = 49.1, right: t = 81.9. 




Figure 5: Same as Figure H] with L = 0.4. The plotting times are now, left: t = 0, middle: t = 8.1, 
right: t = 17.8. 




Figure 6: Same as Figure H] with L = 0.7. The plotting times are now, left: t = 0.07, middle: t = 3.55, 
right: t = 6.77. 
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Figure 7: Different shapes of traveling plateaus. In each figure, the top subplot depicts a zoom 
of the peak, while the bottom subplot displays the whole plateau. Parameter values: L = 0.2, 
a c = a s = D s = 1; a) D c = 0.01, b) D c = 0.2, c) D c = 0.5, d) D c = 0.6. 
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Figure 8: The velocity of the traveling plateau a for different values of L with the parameters a c = 
l,as = 1,D C = 0.1, Ds = 1- We can see that, regardless of the numerical effect, they are almost 
constant. Left: L = 0.1; right: L = 0.4. Top: the full scale, bottom: zoom on the oscillations. 
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